library(Biostrings)
library(XML)

baseDir <- "/home/haseong/alpha/dev/FT-NGS"
targetDir <- "/home/haseong/alpha/dev/FT-NGS"
dbName <- "/home/haseong/alpha/apps/blast/db/nt"
setwd(baseDir)

file.names<-dir(pattern="^KHS.*fna$")
eval.cri<-0.01

for(fn in file.names){
##fn<-file.names[2]
  contig.file<-paste(baseDir, "/", fn, sep="")
  contig.seqs<-readDNAStringSet(contig.file, "fasta")
  
  out.dir<-paste(targetDir, "/out_", fn, sep="")
  contig.file.names<-dir(out.dir, pattern="*.xml")
  bl<-list()
  cnt<-1
  for(cn in contig.file.names){
    score<-c()
    hit.id<-c()
    evalue<-c()
    qfrom<-c()
    qto<-c()
    hfrom<-c()
    hto<-c()
    alen<-c()
    gaps<-c()
    qseq<-c()
    hseq<-c()
    mline<-c()
    hit.obj.name<-c() 
    hit.def<-c()
    
    target.file<-paste(out.dir, "/", cn, sep="")
    doc = xmlTreeParse(target.file, useInternal = TRUE)
    top = xmlRoot(doc)
    ##seq.hits<-top[["BlastOutput_iterations"]][["Iteration"]][["Iteration_hits"]][["Hit"]][["Hit_hsps"]]
    hits<-top[["BlastOutput_iterations"]][["Iteration"]][["Iteration_hits"]]
    idx1<-1
    while(!is.null(hits[[idx1]])){
      sel.hits<-hits[[idx1]]
      hit.hsp<-sel.hits[["Hit_hsps"]]
      idx2<-1
      while(!is.null(hit.hsp[[idx2]])){ 
        hsp<-hit.hsp[[idx2]]
        if(as.numeric(xmlValue(hsp[["Hsp_evalue"]]))<eval.cri){
          hit.obj.name<-c(hit.obj.name, xmlValue(top[["BlastOutput_query-def"]]))
          ##hit.def<-c(hit.def, xmlValue(top[["BlastOutput_iterations"]][["Iteration"]][["Iteration_hits"]][["Hit"]][["Hit_def"]]))
          hit.def<-c(hit.def, xmlValue(sel.hits[["Hit_def"]]))
          hit.id<-c(hit.id, xmlValue(sel.hits[["Hit_id"]]))
          ##hit.len<-xmlValue(seq.hits[[1]][["Hit_len"]])
          score<-c(score, as.numeric(xmlValue(hsp[["Hsp_score"]])))
          evalue<-c(evalue, as.numeric(xmlValue(hsp[["Hsp_evalue"]])))
          qfrom<-c(qfrom, as.numeric(xmlValue(hsp[["Hsp_query-from"]])))
          qto<-c(qto, as.numeric(xmlValue(hsp[["Hsp_query-to"]])))
          hfrom<-c(hfrom, as.numeric(xmlValue(hsp[["Hsp_hit-from"]])))
          hto<-c(hto, as.numeric(xmlValue(hsp[["Hsp_hit-to"]])))
          alen<-c(alen, as.numeric(xmlValue(hsp[["Hsp_align-len"]])))
          gaps<-c(gaps, as.numeric(xmlValue(hsp[["Hsp_gaps"]])))
          qseq<-c(qseq, xmlValue(hsp[["Hsp_qseq"]]))
          hseq<-c(hseq, xmlValue(hsp[["Hsp_hseq"]]))
          mline<-c(mline, xmlValue(hsp[["Hsp_midline"]]))
          idx2<-idx2+1
        }else{
          break;
        }
      }
      idx1<-idx1+1
    }
    if(!is.null(hit.def)){
      bl<-data.frame(cn, hit.def, hit.id, score, evalue, qfrom, qto, hfrom, hto, alen, gaps, qseq, hseq, mline)  
      ##contig.file.names<-dir(out.dir, pattern="*.xml")
      save(bl, file=paste(out.dir, "/", substr(cn, 1, 11), ".RData", sep="")) 
    }
    cat(".");flush.console();
    if(cnt%%50==0){ cat("(", cnt, "/", length(contig.file.names), ") \n");flush.console() }
    cnt<-cnt+1      
  }
  
## end for(fn in file.names){
}




